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Abstract 

We present a new numerical algorithm for evolving the Mixmaster space- 
times. By using symplectic integration techniques to take advantage of the 
exact Taub solution for the scattering between asymptotic Kasner regimes, 
we evolve these spacetimes with higher accuracy using much larger time steps 
than previously possible. The longer Mixmaster evolution thus allowed en- 
ables detailed comparison with the Belinskii, Khalatnikov, Lifshitz (BKL) ap- 
proximate Mixmaster dynamics. In particular, we show that errors between 
the BKL prediction and the measured parameters early in the simulation can 
be eliminated by relaxing the BKL assumptions to yield an improved map. 
The improved map has different predictions for vacuum Bianchi Type IX and 
magnetic Bianchi Type VIo Mixmaster models which are clearly matched in 

the simulation. 
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Mixmaster dynamics (MD), discovered independently by Belinskii, Khalatnikov, and Lif- 
shitz (BKL) ]IJ and Misner 0, describes a system in which pure gravity exhibits chaos (or 
at least a strong sensitivity to initial conditions). Although first discovered in vacuum spa- 
tially homogeneous cosmologies of Bianchi Type IX, MD also occurs in vacuum cosmologies 
of Bianchi Type VIII [JJ] and magnetic Bianchi Type VIo 0. (Other possible arenas for 
MD are given by Jantzen ||.) It is clear that the most general homogeneous cosmology 
exhibits MD and has been conjectured that the same is true in the inhomogeneous case [|TJ. 
Much recent work [g] has focused on the question of whether or not MD is chaotic in any 
invariant sense since computed Lyapunov exponents can be either zero or positive depending 
on the choice of time variable The issue has recently been resolved in favor of chaos 

by Cornish and Levin |H| who have provided a prescription to define the discrete outcomes 
required to exhibit the fractal basins of attraction characteristic of chaos. A parallel issue 



remains, however. BKL (as revised by Chernoff and Barrow JTT| and extended by Berger 
12|| ) derived an approximate MD as a sequence of Kasner models with a map from one Kas- 
ner epoch to the next. So far, the properties of the map have always been consistent with 
those of the full solution to Einstein's equations obtained numerically |7j,P|,|T^] . However, one 
would wish to have more precise criteria for the BKL map's validity and perhaps to measure 
departures from it. In addition, one must untangle the loss of information due to the chaotic 
nature of the dynamics from the accumulation of errors due to finite numerical precision. 

Here, we describe a new numerical algorithm for MD which represents improvement by 
at least two orders of magnitude over standard ODE solvers in the speed with which a 
Mixmaster model may be evolved toward the singularity (without any loss of accuracy). 
We take advantage of the fact that MD as a sequence of Kasners is equivalent to MD as 
a sequence of bounces (transitions between Kasners) and that evolution from one Kasner 
to the next is the exactly solvable Taub cosmology We apply this new method to the 
(diagonal) Bianchi IX and magnetic Bianchi VIo examples of MD in order to compare an 
improved BKL map to the numerical results and to display clearly the role of numerical 
precision in this comparison. 



The models we shall use to illustrate our method are described by the metric 

ds 2 = _ e 2(«H+i) dt 2 + e 2 "^) 2 + e 2 <{a 2 ) 2 + e 2 ^ (a 3 ) 2 . (1) 

Here a, ( and 7 are functions of t, and a\, a 2 and 03 are time independent, orthogonal 
forms, invariant under the group of spatial symmetries. The dynamics of these spacetimes 
is therefore just the determination of the logarithmic scale factors (LSFs) a, ( and 7 as 
functions of t. The Einstein equations for the models may be obtained by variation of the 
Hamiltonian H = Hk + H p (and the constraint H = 0) where 

2 H k = 3 (p 2 a + p\ + p 2 ) - 6 {p a p ( + p a p^ + p c p 7 ) . (2) 

The potential term H p is a function of the LSFs and depends on the type of homogeneous 
cosmology being treated. For vacuum Bianchi type IX or magnetic Bianchi Type VIo, we 
have 

Hp = c 2 e 2ba + e 4? + e 4 ^ - 2 (ae 2 ^ + ae 2 ^' + de 2 ^) . (3) 

Here a = 1, 6 = 2, c — 1 and d = 1 for vacuum Bianchi Type IX, while a = 0, b = 1, 
c = a/£ and d — — 1 for magnetic Bianchi Type VIo, and £ is a constant that depends on 
the strength of the magnetic field. Note that the magnetic field in Type VIo provides the 
potential wall that is furnished by curvature in Type IX. When the logarithmic scale factors 
are large and negative, the potential term H p is negligible so that the dynamical system is 
then essentially a free particle: a, ( and 7 are linear functions of t. This is a Kasner epoch. 
The bounces between epochs occur in those periods when the potential is not negligible. 
Since the potential consists of terms that are exponentials of the LSFs, each of the bounces 
occurs in a brief time period (in terms of t). 

To evolve the Mixmaster spacetimes numerically we use the method of symplectic integra- 



tion |TjJT^]. This method begins with a dynamical system with Hamiltonian H = Hi + H 2 
where the equations obtained by variation of Hi and H 2 are separately exactly solvable. 
For a Hamiltonian H let U(H, At) be the operator that evolves the system for a time At. 
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Consider the operator 

U 2 (H, At) = U{H U At/2)U(H 2 ,At) U(H ± , At/2) . (4) 

then to second order in At the operator U 2 (H, At) agrees with U (H, At). Using the operator 
U 2 (H,At) one can, by iteration, construct evolution operators that agree with U(H,At) to 
any desired order. Explicitly the 2n + 2 order approximation to U(H, At) is 

U 2n+2 (H, At) = U 2n (H, s n At) U 2n (H, (1 - 2s n )At) U 2n (H, s n At) (5) 

where s n = 1/(2 - 2 1 /( 2n+1 )). Note that though evolution using U 2n is formally 2n order 
accurate, the error may be much smaller than one would expect. Recall that if H 2 = then 
U 2 (and therefore U 2n ) is exactly equal to U. Similarly if the evolution takes place in a region 
of phase space where H 2 is extremely small then U 2n is an extremely good approximation 
to U, a much better approximation than, e.g. 2n order Runge-Kutta. Thus in such regions 
of phase space one can take very large time steps without introducing large inaccuracy. 

To apply the symplectic integration method to the homogeneous cosmologies it is neces- 
sary to divide their Hamiltonian into two pieces, each of which is exactly solvable. It turns 
out that H k and H p are each exactly solvable since H k is a free particle Hamiltonian that 
yields the Kasner solution while H p contains no momenta and thus yields trivial equations 
of motion. However, splitting the Hamiltonian in this way still necessitates the use of ex- 
tremely small time steps at the bounces. To take full advantage of the symplectic method 
we would like to split the Hamiltonian in such a way that H 2 is very small in the region 
of phase space in which the evolution is taking place. Note that the largest contribution 
to the potential term is always of the form c 2 e 2ba where b and c are the constants defined 
in (3) for a corresponding to a magnetic wall and c = 1, b — 2 for a corresponding to a 
curvature wall. Call the other two LSFs ( and 7. Denote by p a , and p 7 the momenta 
corresponding to a, ( and 7 respectively. We split the Hamiltonian into Hi = H k + c 2 e 2ba 
and H 2 = H, p — c 2 e 2ba . But Hi is the Hamiltonian for the exactly solvable Taub model while 
H 2 still contains no momenta. At each time step, the code identifies the largest LSF and 
implements the appropriate split. 
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In our notation, the Taub solution is the following: Since Hi is independent of ( and 7 
it follows that p^ and p 7 are constants. The remaining variables evolve as follows: Define 
q, k and r by q = — a(t), 

k =[q 2 + 6c 2 e 2b ^] 1/2 , (6) 

r = In (^cosh kbAt + | sinh £;6At^ . (7) 

Then we have 

a(t + At) = a(t) - T - , (8) 

p a (t + At) = p Q (t) - — e 265(t) ~ r sinhA;6At , (9) 

C(t + At) = C(t) + £ - Qp.At , (10) 

7(t + At) =7(t) + T - - 6p c At . (11) 

Note that these formulas give exact solutions of the dynamics of ifi. There is no approxi- 
mation and no assumption that At is "small." The dynamics of H 2 are trivial. Since H 2 is 
independent of p a , and p 7 it follows that a, ( and 7 are constants. The evolution of p a 
is then given by 

p a (t + At)=p Q (t) - ^ At (12) 

and correspondingly for the other momenta. 

Our computer code implements the sixth order symplectic integration algorithm with 
this split of the Hamiltonian. The time step is increased if evolution with At and At/2 
yields the same solution to a desired accuracy. For comparison we also evolved the same 
Mixmaster spacetimes with a fourth order, adaptive stepsize Runge-Kutta code and a sixth 
order symplectic algorithm with H = H k + H p . A comparison of the results of the three 
codes for type IX is shown in Fig. 1 for the region around a typical bounce. Note that the 
symplectic integration code achieves accurate evolution in approximately one hundred times 
fewer steps than the Runge-Kutta code. The key here is that any adaptive step size algorithm 
in any ODE solver can accurately reproduce the Kasner solution. However, in standard 
methods, At must be drastically decreased when a potential term becomes non- negligible. 



If care is not taken, a too large time step will cause overshoot of the potential wall leading 
to numerical disaster. In the new algorithm, the largest term in (3) is identified and the 
appropriate Taub solution used to effortlessly propagate the trajectory through the bounce. 
There is no need to decrease the time step. One can easily reach \Q\ = ||a + C + 7| ~ 10 35 in 
only 50,000 time steps corresponding to more than 150 epochs. Previous simulations have 
reported results for fewer than 30 epochs with < 10 8 while requiring significantly more 
computer time |§. If the constraint is not enforced at all, none of the algorithms discussed 
here can evolve more than about 15 epochs with the initial data used here. To maintain 
accuracy, it is necessary to enforce the Hamiltonian constraint at least every few time steps 



depending on the precision required. Recently, Gundlach and Pullin [[TJ| have argued that 
the constraints must be enforced in generic numerical relativity codes. Our new algorithm 
allows simulations to be run for a sufficiently large number of epochs that the need to enforce 
the constraint here is also seen. 

Each Kasner epoch can be specified by the values of four parameters: (u,v,pn, k) . As 
shown by several authors [|IT|,|Tf these parameters can be defined throughout the evolution 
of a Mixmaster spacetime, though they are (approximately) constant only within a Kasner 



epoch. Following reference JH]] we define pn = p a + p^ + p 1 , and the Kasner indices p 1 = 
-d/(3pn) , P2 = ~C/(3pn) and p 3 = -7/(3pn)- Now reorder the pi so that pi < p 2 < Pz- 
Let a, ( and 7 be respectively the LSFs corresponding to pi, p 2 and p 3 . Define u, v and k 
by 

u = -l - ?* , (13) 
Pi 

P3 (pxC - p 2 a, 



k = P3 (£---) ■ (15) 

\P2 Pi J 

BKL made the following approximation to MD: (i) between bounces the spacetime is exactly 
Kasner, (ii) the bounce occurs instantaneously and (iii) at a bounce one of the LSFs vanishes. 
This approximation yields a rule (the BKL map) that gives the values of the parameters 
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in the n + 1st epoch in terms of their values in the nth epoch. It is well-known || that 
these assumptions are less valid early in the simulation and improve as the singularity is 
approached. For pq, the BKL rule is 



/ u n - u n + 1 > 

Pn,n+1 =Pn,n -5— ; — • (16) 



For the other parameters, the BKL map depends on whether u n is greater than 2. For u n > 2, 
the BKL map is u n+ i — u n — 1, fn+i = t> n + 1 and K n+ i = k„. For w n < 2 (referred to as 
an era change), the rule is u n+ i = 1/ (u n - 1), v n +i = 1 + (l/v n ) and = v n K n / (u n - 1). 
However, the known exact dynamics of H\ provides a better approximation to MD than 
does this BKL approximate dynamics — i.e. the sequence of bounces is a better description 
than the sequence of Kasners. It therefore yields a better map of the parameters (u, v,pq,k) 
from one epoch to the next. This "improved BKL map" actually agrees with the BKL map 
for u and pn but adds corrections to v and k. Let the quantity w be given by 

2 , / c u 2 + u + l\ 

w = T In -= 17 

b \V6 up n J 

where b and c are defined as in (3) for the vacuum Type IX and magnetic VIo cases. Within 
any era the improved BKL map for v and k is 

V n+l = V n + 1 + — — r , 18 

1 + (Kn/Wn) 

K n+1 = K n + W n , (19) 

while for an era change the improved BKL map gives 

V n + [U n + 1) {Wn/l^n) 

V n K n + (U n + l)w n 

K n +1 = : • (21) 

U n - 1 

Note that in the limit of vanishing w n / K n the improved BKL map agrees with the BKL 
map. As the singularity is approached, \n n \ —>■ oo so that the BKL terms overwhelm the 
corrections. In Fig. 2, 20 epochs of a typical magnetic Bianchi VIq trajectory are displayed in 



the anisotropy plane. The axes are |L8| @±/\Q\ where (3+ = |(C + 7 — 2a), /3_ = ^/§(C + 7) so 
that the exponential potential terms have fixed locations. We have chosen the vertical wall 
to be due to the magnetic field (£e 2a ) with the other two walls (e 4( % e 4 " 7 ) due to curvature. 
The epoch numbers are identical to those in Fig. 3 which shows the fractional difference 
between measured and predicted (on the basis of the previous measured value) values of v 
for the BKL and improved maps for bounces off magnetic (VI) and curvature (IX) walls. It 
is clearly seen that extremely accurate agreement is obtained alternately for the magnetic 
and curvature improved maps as the trajectory correspondingly bounces off the respective 
walls. A new era with bounces between two curvature walls is also seen in epochs 12-14. 
Fig. 4 shows a later stage in the evolution when all three maps converge as expected. The 
results for the parameter k are identical. 

Finally, we consider tracking an initial value of u through almost 200 epochs with the 
Hamiltonian constraint enforced at every time step. The sequence is begun at a point where, 
during an era, the ODE solution follows the BKL map for u to the desired precision. Unlike 
the predictions previously described, loss of information and numerical error are allowed 
to accumulate. Mathematica [T9| is used to obtain the sequence of BKL w-values to 60 



significant figures (SFs) while our ODE solver is run in quadruple precision (32 SFs). Integer 
parts for u in the latter different from those predicted by the former arise in our simulation 
at epoch 161, precisely where a 32 SF Mathematica sequence completely loses precision due 
to era changes. In fact, examination of the u sequence shows that epoch 161 begins the 32nd 
era, as expected. (Similarly, the double precision simulation deviates at epoch 47, the start 
of the 16th era.) Given initial data specified to N SFs, the numerical simulation follows the 
true Mixmaster trajectory through « N eras, until the initial information is lost. However, 
as can be deduced from Fig. 4, the numerical solution remains close to some true Mixmaster 
trajectory for any sequence of N eras in the simulation. We emphasize here that these effects 
are not visible early in the simulation. Our new algorithm allows easy study of Mixmaster 
models for more than one hundred epochs to permit careful examination of issues relating 
to the validity of the (improved) BKL approximation. 



ACKNOWLEDGMENTS 



This work was supported in part by National Science Foundation Grants PHY9507313 
(BKB) and PHY9408439 (DG) to Oakland University (OU). DG was also supported in part 
by a Cottrell College Science Award of Research Corporation to OU. BKB would like to 
thank the Institute for Geophysics and Planetary Physics at Lawrence Livermore National 
Laboratory and the Department of Astronomy at the University of Michigan for hospitality. 
Part of this work was performed by ES in partial fulfillment of the requirements for the 
M.S. in Physics degree at OU. 



9 



[1] 

[2] 
[3] 
[4] 
[5] 
[6] 



REFERENCES 

Belinskii V A, Lifschitz E M, and Khalatnikov I M 1971 Usp. Fiz. Nauk. 13 463 
Misner C W 1969 Phys. Rev. Lett 22 1071 
Halpern P 1987 J. Gen. Rel. Grav. 19 73 

LeBlanc V G , Kerr D, and Wainwright J 1995 Class. Quantum Grav. 12 513 
Jantzen R T 1986 Phys. Rev. D 33 2121 

See articles by Hobill D W, Misner C W, Berger B K, Burd A, Francisco G and Matsas 
G E A, Rugh S E, and Ma P K-H and Wainright J in Hobill D, Burd A, and Coley, A 
(eds) 1994 Deterministic Chaos in General Relativity (New York: Plenum) 

Rugh S E 1990 Cand. Scient. Thesis, Niels Bohr Inst.; Rugh S E, Jones B J T 1990 
Phys. Lett. A147 353 

Berger B K 1991 Gen. Rel. Grav. 23 1385 

Ferraz K, Francisco G, Matsas G E A 1991 Phys. Lett. 156A 407 
Cornish N J, Levin J J kr-qc/9605029 



[9 
[10 

[11 
[12 

[13 

[14 

[15 

[16 

[17] Berger B K 1994 Phys. Rev. D 49 1120 

[18] Moser A R, Matzner R A, Ryan M P Jr 1973 Ann. Phys. (N. Y.) 79 558 

10 



ChernofT D F, Barrow J D 1983 Phys. Rev. Lett. 50 134 
Berger B K 1993 Phys. Rev. D 47 3222 
Taub A 1951 Ann. Math. 53 472 

Fleck J A, Morris J R, and Feit M D 1976 Appl. Phys. 10 129 
Suzuki M 1990 Phys. Lett. A146 319 
Gundlach C, Pullin J |gr-qc/ 9606022 



[19] Mathematica, Wolfram Research 



FIGURES 

FIG. 1. Comparison of algorithms for MD. A typical Mixmaster bounce is shown in the 
anisotropy plane. Crosses indicates every 10th point on a 4th order Runge-Kutta evaluation of the 
trajectory, while circles indicate every 10th point for a 6th order standard symplectic evaluation. 
The filled squares indicate every point using the new algorithm. The inset shows the details closest 
to the bounce. 

FIG. 2. The first 20 epochs of a typical magnetic Bianchi VIo trajectory in the rescaled 
anisotropy plane where the potential is shown as a fixed equilateral triangle. Some epoch numbers 
are shown. The vertical potential wall is produced by the magnetic field while the others are due 
to curvature. 

FIG. 3. Fractional difference between measured and predicted values of v vs. epoch number 
using the BKL prediction (crosses), the curvature wall prediction (squares), and the magnetic wall 
prediction (circles). Triangles mark the first epoch of an era. 

FIG. 4. Convergence of the BKL and improved maps as the evolution shown in Figures 2 and 
3 is continued toward the singularity. The symbols are the same as those in Fig. 3. 
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